A Monotonically Convergent Newton Iteration for the Quantiles of any Unimodal Distribution, with Application to the Inverse Gaussian Distribution

One of the most basic and commonly used numerical computations in probability and statistics is to evaluate the random deviate corresponding to any given tail probability for a given probability distribution. The deviate corresponding to a given probability is called the quantile. Quantiles usually need to be computed by numerical approximation, and the need often arises to compute quantiles for probability distributions for which reliable code is not readily available. The purpose of this article is to point out a simple but elegant result that applies to all continuous unimodal distributions. Newton’s method for finding the quantiles of a continuous unimodal distribution is always monotonically convergent when started from the mode of the distribution. This provides a simple, accurate and numerically reliable method of computing quantiles for any continuous unimodal distribution, given that the cumulative distribution and probability density functions can be evaluated accurately. The monotonic Newton iteration has been implemented in the qinvgauss function of the R package statmod to compute quantiles of inverse Gaussian distributions. The resulting function proves to be faster, more accurate and more reliable than existing functions for the same purpose, even without sophisticated optimization.